Volatility-Surface-without-Interpolation.R

# Required libraries
library(stats)
library(pracma)

# Black-Scholes call price formula
black_scholes_call <- function(S, K, T, r, sigma) {
  d1 <- (log(S / K) + (r + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
  d2 <- d1 - sigma * sqrt(T)
  call_price <- S * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
  return(call_price)
}

# Function to find the implied volatility
implied_volatility <- function(S, K, T, r, market_price) {
  objective_function <- function(sigma) {
    black_scholes_call(S, K, T, r, sigma) - market_price
  }
  
  # Use uniroot to find the implied volatility
  result <- tryCatch({
    uniroot(objective_function, lower = 0.001, upper = 5)$root
  }, error = function(e) NA)
  
  return(result)
}

# Given data
S <- 100  # Current index price
r <- 0.03  # Risk-free interest rate

# Market prices matrix
market_prices <- matrix(c(
  40.2844, 42.4249, 50.8521, 59.1664,
  30.5281, 33.5355, 42.6656, 51.2181,
  21.0415, 24.9642, 34.4358, 42.9436,
  12.2459, 16.9652, 26.4453, 34.789,
  5.2025, 10.1717, 19.4706, 27.8938,
  1.3448, 5.4318, 14.4225, 23.3305,
  0.2052, 2.7647, 11.2103, 20.7206,
  0.0216, 1.4204, 9.1497, 19.1828,
  0.0019, 0.7542, 7.741, 18.1858
), nrow = 9, byrow = TRUE)

# Strikes and maturities
strikes <- c(60, 70, 80, 90, 100, 110, 120, 130, 140)
maturities <- c(0.25, 0.5, 1, 1.5)

# Create a matrix to store implied volatilities
implied_vols <- matrix(NA, nrow = length(strikes), ncol = length(maturities),
                       dimnames = list(paste0("Strike_", strikes), paste0("Maturity_", maturities)))

# Calculate implied volatilities for all combinations
for (i in seq_along(strikes)) {
  for (j in seq_along(maturities)) {
    implied_vols[i, j] <- implied_volatility(S, strikes[i], maturities[j], r, market_prices[i, j])
  }
}

# Print the implied volatilities matrix
print(implied_vols)

##            Maturity_0.25 Maturity_0.5 Maturity_1 Maturity_1.5
## Strike_60             NA    0.5865828  0.8175021    0.9478482
## Strike_70      0.2486434    0.5220188  0.7158389    0.8294040
## Strike_80      0.3148182    0.4545265  0.6174672    0.7155847
## Strike_90      0.2801724    0.3894853  0.5277066    0.6147351
## Strike_100     0.2426795    0.3368458  0.4608755    0.5465596
## Strike_110     0.2173596    0.3070502  0.4299974    0.5243019
## Strike_120     0.2056962    0.2975798  0.4272046    0.5335653
## Strike_130     0.2022396    0.2989333  0.4376976    0.5558278
## Strike_140     0.2030054    0.3049178  0.4531255    0.5819529

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# Prepare the data for plotly
df <- expand.grid(Strike = strikes, Maturity = maturities)
df$Volatility <- as.vector(implied_vols)

plot_ly(data = df, x = ~Maturity, y = ~Strike, z = ~Volatility, type = "mesh3d") %>%
  layout(scene = list(
    xaxis = list(title = "Maturity"),
    yaxis = list(title = "Strike"),
    zaxis = list(title = "Implied Volatility")
  ))